(* Mathematica package *)


(*Internal implementation, less restrictive and more prone to input error.  
To be used internally with parent functions passing correct arguments*)
AbelesS[pol_, omega_, thetalist_, nlist_, Dlist_] := 
	Dot @@ (Riffle[
		RefractionMv[pol, Most[thetalist], Rest[thetalist], Most[nlist], Rest[nlist]],
		PhaseMv[omega, thetalist[[2 ;; -2]], nlist[[2 ;; -2]], Dlist]]);
      
(*End User definition of S, more restrictive pattern-match for args, executes ~ x1.09 as internal definition. *)

(*Compute Refraction Angles Internally*)
AbelesS[pol_?PolQ, omega_, Theta0_, n0_, NvDv:{{_,_}..}, nkp1_] :=
	Block[{ThetaLIST = {Theta0} ~Join~ RefractionAngle[Theta0, n0, NvDv[[All,1]]~Join~{nkp1}],
	       nLIST = {n0} ~Join~ NvDv[[All,1]] ~Join~ {nkp1}},
		Dot @@ (Riffle[
			RefractionMv[pol, ThetaLIST[[1 ;; -2]], ThetaLIST[[2 ;;]],nLIST[[1 ;; -2]], nLIST[[2 ;;]]],
			PhaseMv[omega, ThetaLIST[[2;;-2]], Sequence @@ Transpose[NvDv]]])];

(*Precompute refraction angles*)
AbelesS[pol_?PolQ, omega_, ThetavNvDv : {{_, _}, {_, _, _} .., {_, _}}] :=
	Dot @@ (Riffle[
		RefractionMv[pol, ThetavNvDv[[1 ;; -2, 1]], ThetavNvDv[[2 ;;, 1]],ThetavNvDv[[1 ;; -2, 2]], ThetavNvDv[[2 ;;, 2]]],
		PhaseMv[omega, Sequence @@ Transpose[ThetavNvDv[[2 ;; -2]]]]]);

AbelesS[___]/;Message[General::badargs,AbelesS] :=
	"dummy result, so symbol returned unevaluated";